Dado un sistema masa resorte amortiguador como el de la figura siguiente:
Graficar la trayectoria del sistema, y animar los componentes fisicos de manera esquematica.
Recordando la práctica anterior, podemos obtener el comportamiento del sistema:
$$ F(t) - k x(t) - c \dot{x}(t) = m \ddot{x}(t) $$a traves de su función de transferencia:
$$ G(s) = \frac{\frac{1}{m}}{s^2 + \frac{c}{m} s + \frac{k}{m}} $$
In [1]:
from control import step, tf
In [2]:
m = 1200
k = 15000
c = 1500
F = 1
In [3]:
G = tf([0, 0, 1/m], [1, c/m, k/m])
G
Out[3]:
In [4]:
y, t = step(G)
Una vez que tenemos los datos, nos disponemos a graficarlos, ahora utilizaremos un estilo de matplotlib, para ver los estilos disponibles, tan solo tenemos que preguntar:
In [5]:
%matplotlib inline
from matplotlib.pyplot import figure, style, plot
style.use("ggplot")
In [6]:
style.available
Out[6]:
Ahora solo damos el tamaño de la imagen, y obtenemos el objeto del axis de la grafica con el metodo gca para poder graficar:
In [7]:
f = figure(figsize=(8, 8))
axi = f.gca()
axi.plot(t, y);
Sin embargo el comportamiento del sistema a penas es visible, $0.00006 m = .6 mm$ en posición es practicamente nada. Recordemos que usamos tan solo usamos una fuerza de $1 N$, para amplificar esta entrada del sistema utilizamos otro bloque con una ganancia de la magnitud de la fuerza a aplicar:
In [8]:
K = tf([1000], [1])
K
Out[8]:
In [9]:
y, t = step(K*G)
f = figure(figsize=(8, 8))
axi = f.gca()
axi.plot(t, y);
Muy bien, nuestra suspensión se ha movido $10 cm$ con una fuerza comparable con el peso de una persona. Ahora notamos que a la derecha tenemos espacio vacio, por lo que pondremos limite al eje $x$:
In [10]:
f = figure(figsize=(16, 8))
axi = f.gca()
axi.plot(t, y)
axi.set_xlim(0, 10);
Y la suspensión no se mueve de abajo hacia arriba, sino al reves, multiplicamos por $-1$ los datos de $y$:
In [11]:
f = figure(figsize=(16, 8))
axi = f.gca()
axi.plot(t, -y)
axi.set_xlim(0, 10);
Para poder ver en una escala razonable esta gráfica pondremos un limite al eje $y$ tambien:
In [12]:
from numpy import linspace
In [13]:
t = linspace(0, 10, 1000)
y, t = step(K*G, t)
f = figure(figsize=(16, 8))
axi = f.gca()
axi.plot(t, -y)
axi.set_xlim(0, 10)
axi.set_ylim(-1.2, 0);
Si ahora metemos todo este codigo dentro de una función para graficar, dependiendo de los datos del problema, se verá asi:
In [14]:
def suspension(m, c, k, F):
G = tf([0, 0, 1/m], [1, c/m, k/m])
K = tf([F], [1])
t = linspace(0, 10, 1000)
y, t = step(K*G, t)
f = figure(figsize=(16, 8))
axi = f.gca()
axi.plot(t, -y)
axi.set_xlim(0, 10)
axi.set_ylim(-1.2, 0);
# La siguiente linea devuelve los resultados numericos de la simulacion,
# si pones un punto y coma (;) despues de usar la funcion, o si guardas
# en variables, podras ver la grafica solamente
return t, y
In [15]:
suspension(1200, 1500, 15000, 100*9.81);
¿Cuanto crees que pese una bola de demolición? ¿$500 kg$? ¿$1000 kg$? Grafica el comportamiento del sistema ante una fuerza de estas magnitudes, utilizando la función que acabamos de crear.
In [29]:
t1, x1 = suspension(1200, 1500, 15000, ) # Agrega la fuerza necesaria para la primer simulacion
In [31]:
t2, x2 = # Completa el codigo para la segunda simulacion
In [32]:
from pruebas_3 import prueba_3_1
prueba_3_1(t1, x1, t2, x2)
Esta es la parte divertida, para interactuar con los datos del problema, podemos usar una función especifica de IPython, tan solo tenemos que dar una función, y los rangos de variación de sus argumentos, para tener:
In [ ]:
# Se importan widgets de IPython para interactuar con la funcion
from IPython.html.widgets import interact, fixed
In [ ]:
# Se llama a la funcion interactiva
interact(suspension, m=(1000, 2500), c=(0, 2500), k=(0, 20000), F=(1, 10000));
A quien no le gustan los GIF's?
Hasta Chuck Norris lo sabe...
En fin, hoy haremos nuestro GIF... del comportamiento del sistema!
Pero antes, tenemos que definir como graficaremos nuestros elementos, empecemos con el resorte. Basicamente en una linea retorcida, por lo que podemos simplemente dar las coordenadas y decirle a plot que haga una linea entre estas:
In [ ]:
from numpy import zeros, arange
In [ ]:
coordenadas_resorte = zeros(20)
In [ ]:
coordenadas_resorte[5:15] = 0.5*(-1)**arange(10)
In [ ]:
coordenadas_resorte
In [ ]:
plot(coordenadas_resorte, arange(20), lw = 2, color="gray")
Ya que tenemos esto, podemos meterlo en una funcion que nos de las coordenadas y el codigo nos queda mas simple:
In [ ]:
def coordenadas_resorte(y0, yf, x, N):
ancho = (y0 - yf)/5
xs = zeros(N) + x
xs[N/4:N*3/4] = ancho*0.5*(-1)**arange(N/2) + x
return xs, linspace(y0, yf, N)
In [ ]:
x, y = coordenadas_resorte(-5, 5, 0, 20)
f = figure(figsize=(8, 8))
axi = f.gca()
axi.plot(x, y, lw=2, color="gray")
axi.set_xlim(-10, 10)
axi.set_ylim(-10, 10);
El amortiguador es un poco mas complejo, son 3 lineas, 2 para las patitas y una para el cuerpo:
In [ ]:
y0 = -5
yf = 5
alto = yf - y0
cuerpo = alto/2
pata = alto/4
ancho = (yf - y0)/5
x = 2
f = figure(figsize=(8, 8))
axi = f.gca()
axi.plot([-0.5*ancho + x, -0.5*ancho + x, 0.5*ancho + x, 0.5*ancho + x],
[yf - pata, y0 + pata, y0 + pata, yf - pata],
[x, x], [y0 + cuerpo, yf],
[x, x], [y0, y0 + pata], lw=2, color="gray")
axi.set_xlim(-10, 10)
axi.set_ylim(-10, 10);
De la misma manera lo metemos en una función, para que nuestro codigo quede:
In [ ]:
def coordenadas_amortiguador(y0, yf, x):
ancho = (y0 - yf)/5
alto = yf - y0
cuerpo = alto/2
pata = alto/4
ax = [-0.5*ancho + x, -0.5*ancho + x, 0.5*ancho + x, 0.5*ancho + x]
ay = [yf - pata, y0 + pata, y0 + pata, yf - pata]
bx = [x, x]
by = [y0 + cuerpo, yf]
cx = [x, x]
cy = [y0, y0 + pata]
return ax, ay, bx, by, cx, cy
In [ ]:
ax, ay, bx, by, cx, cy = coordenadas_amortiguador(-4, 6, 1)
f = figure(figsize=(8, 8))
axi = f.gca()
axi.plot(ax, ay, bx, by, cx, cy, lw=2, color="gray")
axi.set_xlim(-10, 10)
axi.set_ylim(-10, 10);
Por ultimo, para la masa tan solo necesitamos un rectangulo, el cual tenemos que importar de la libreria de matplotlib del modulo patches:
In [ ]:
from matplotlib.patches import Rectangle
In [ ]:
f = figure(figsize=(8, 8))
axi = f.add_subplot(111, autoscale_on=False, xlim=(-3, 3), ylim=(-3, 3))
carro = Rectangle((0,0), 1, 1, lw=0.5)
carro.set_xy((-0.5, -0.5))
axi.add_patch(carro);
Y uniendo todo, nos queda el codigo:
In [ ]:
f = figure(figsize=(8, 8))
axi = f.add_subplot(111, autoscale_on=False, xlim=(-2, 2), ylim=(-3, 1))
carro = Rectangle((0,0), 1, 1, lw=0.5)
carro.set_xy((-0.5, -3))
axi.add_patch(carro)
x, y = coordenadas_resorte(0, -2, -0.25, 20)
axi.plot(x, y, lw=2, color="gray")
ax, ay, bx, by, cx, cy = coordenadas_amortiguador(-2, 0, 0.25)
axi.plot(ax, ay, bx, by, cx, cy, lw=2, color="gray");
axi.set_ylim(-4, 0);
Por ultimo, para animar esto, tenemos que obtener los puntos de la trayectoria:
In [ ]:
m = 1200
c = 1500
k = 15000
F = 4000
G = tf([0, 0, 1/m], [1, c/m, k/m])
K = tf([F], [1])
t = linspace(0, 10, 1000)
y, t = step(K*G, t)
Importamos el modulo de animacion de la libreria matplotlib:
In [ ]:
from matplotlib import animation
Y utilizamos el siguiente codigo para la creacion de cada uno de los cuadros del GIF:
In [ ]:
# Se define el tamaño de la figura
fig = figure(figsize=(8, 8))
# Se define una sola grafica en la figura y se dan los limites de los ejes x y y
axi = fig.add_subplot(111, autoscale_on=False, xlim=(-1.5, 1.5), ylim=(-2, 1))
# Se hacen invisibles los ejes
axi.axes.get_xaxis().set_visible(False)
axi.axes.get_yaxis().set_visible(False)
# Se quitan los bordes de la grafica
axi.axes.spines["right"].set_color("none")
axi.axes.spines["left"].set_color("none")
axi.axes.spines["top"].set_color("none")
axi.axes.spines["bottom"].set_color("none")
# Se utilizan graficas de linea para el resorte y amortiguador
resorte, = axi.plot([], [], lw=2, color='gray')
amor1, = axi.plot([], [], lw=2, color='gray')
amor2, = axi.plot([], [], lw=2, color='gray')
amor3, = axi.plot([], [], lw=2, color='gray')
# Se utiliza un rectangulo para la masa
masa = Rectangle((10,10), 1, 1, lw=0.5)
def init():
# Esta funcion se ejecuta una sola vez y sirve para inicializar el sistema
resorte.set_data([], [])
amor1.set_data([], [])
amor2.set_data([], [])
amor3.set_data([], [])
masa.set_xy((0, 0))
axi.add_patch(masa)
return resorte, amor1, amor2, amor3, masa
def animate(i):
# Esta funcion se ejecuta para cada cuadro del GIF
# Se obtienen las coordenadas del resorte y se meten los datos en su grafica de linea
xs, ys = coordenadas_resorte(-y[i], 1, -0.25, 20)
resorte.set_data(xs, ys)
# Se obtienen las coordenadas del amortiguador y se meten los datos en cada una de
# las graficas de linea
ax, ay, bx, by, cx, cy = coordenadas_amortiguador(-y[i], 1, 0.25)
amor1.set_data(ax, ay)
amor2.set_data(bx, by)
amor3.set_data(cx, cy)
# Se meten los datos de ubicacion de la masa
masa.set_xy((-0.5, -1 - y[i]))
return resorte, amor1, amor2, amor3, masa
# Se hace la animacion dandole el nombre de la figura definida al principio, la funcion que
# se debe ejecutar para cada cuadro, el numero de cuadros que se debe de hacer, el periodo
# de cada cuadro y la funcion inicial
ani = animation.FuncAnimation(fig, animate, arange(1, len(y)), interval=25,
blit=True, init_func=init)
# Se guarda el GIF en el archivo indicado
ani.save('./imagenes/masa-resorte-amortiguador.gif', writer='imagemagick');